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Abstract. We review the current status of attempts to numerically model the merger of 
neutron star-neutron star (NSNS) and black hole-neutron star (BHNS) binary systems, and 
we describe the understanding of such events that is emerging from these calculations. To 
accurately model the physics of NSNS and BHNS mergers is a difficult task. It requires 
solving Einstein's equations for dynamic spacetimes containing black holes. It also requires 
evolving the hot, supernuclear-density neutron star matter together with the magnetic and 
radiation fields that can influence the post-merger dynamics. Older studies concentrated on 
either one or the other of these challenges, but now efforts are being made to model both 
relativity and microphysics accurately together. These NSNS and BHNS simulations are then 
used to characterize the gravitational wave signals of such events and to address their potential 
for generating short-duration gamma ray bursts. 



PACS numbers: 04.25D-, 04.25.dk, 04.30.Db, 47.75.+f, 95.30.Sf 

Compact object binaries consisting of either two neutron stars (NSNS binaries) or a 
black hole (BH) and a neutron star (NS) (BHNS binaries) are unique laboratories for studying 
both strongly-curved spacetime and supernuclear-density matter. Such binaries radiate away 
orbital energy through gravitational waves, causing the two objects to slowly spiral in closer 
to one another and, eventually, to collide and merge. 

Estimating the rate of NSNS/BHNS mergers is not easy. In our galaxy, there are several 
observed NSNS binaries that should merge due to gravitational radiation-induced inspiral in 
less than a Hubble time. (The double pulsar system J0737-3039 will coalesce in 85Myr HI.) 
Extrapolating from these observations, the rate of NSNS mergers per Milky Way equivalent 
galaxy has been estimated in the range 40-700 per Myr [2J. In contrast, no BHNS binaries 
have yet been detected. The merger rate can also be estimated theoretically. NSNS and BHNS 
binaries are thought to form from binary systems of two massive (> 8M ) stars. Population 
synthesis calculations model the evolution of a large number of massive stellar binaries to 
determine the likelihood of producing tight NSNS or BHNS systems (see e.g. [J3j SI for 
details.). There are large uncertainties regarding several aspects of the evolution of these 
binaries, leading to more than an order of magnitude uncertainty in the predicted NSNS and 
BHNS merger rates. Nevertheless, the predicted NSNS merger rate is consistent with that 
inferred from the number of observed systems, and the BHNS rate is predicted to be about 
two orders of magnitude lower Q. 
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NSNS/BHNS mergers might have a number of observable effects. The most direct way 
to observe them is through the gravitational waves emitted during the late inspiral and merger. 
Population synthesis calculations predict that Advanced LIGO should observe of order 10 
NSNS merger events per year, and about one BHNS event 0. These events are also thought 
to be the most promising explanation of the majority of observed short-duration gamma ray 
bursts (GPvBs). Most GRB models require a BH surrounded by a large accretion disk. Energy 
is extracted, either from the disk or the black hole's spin, and is dumped into a baryon- 
poor region around the poles. This drives the ultrarelativistic outflow needed to explain 
the observed athermal spectra of GRBs. (For a review of short GRB progenitor models, 
see [61.) Only numerical simulations can test the assumptions of this model: a massive post- 
merger accretion disk, a baryon-poor region, and efficient energy extraction through neutrinos 
or magnetic fields. Finally, NSNS/BHNS mergers may be important for understanding the 
observed abundances of the heavy elements that are formed by rapid neutron capture in the 
r-process 0. (See (8]| for a recent review of the r-process and its proposed sites.) Here again 
simulations are needed to determine how much and what kind of matter is ejected during such 
mergers. 

In the next section, we review the physics and numerics involved in modeling 
NSNS/BHNS mergers. Then, for each type of binary, we describe the current understanding 
of the merger process derived from simulations. For another useful review of this field, see 
Faber 0. 

1. The numerical challenge 

Neutron stars have compactions in the range GM/Rc 2 ~ 0.1 — 0.2, so we expect general 
relativistic (GR) effects to be important to their evolution. Important studies have nevertheless 
been carried out using Newtonian physics (see below). For BHNS binaries, this means the 
BH is modeled as a Newtonian point mass. A next step towards GR is to modify the point 
mass potential to mimic certain GR effects, such as having an innermost stable circular orbit 
(ISCO) for test particle orbits ifTOl [TT| . In many cases, a good approximation to full GR 
is the conformal flatness approximation (CFA) lfT2l . CFA assumes the spatial metric g,ij is 
conformally flat at all times: = ip^Sij. (Throughout this paper, Latin indices run from 1 to 
3, while Greek indices run from to 3, with being the time index.) From this assumption, 
one can derive elliptic equations for tp (from the Hamiltonian constraint), for the shift ft 1 (from 
the condition that remain conformally flat together with the momentum constraint), and for 
the lapse a (from the maximal slicing condition). CFA is always accurate to first PN order; 
for a spherically- symmetric system it is exact, because for such systems, there are always 
coordinates in which g^ is conformally flat. For spinning NS, errors from assuming CFA are 
a few percent lfT3~l . The spatial metric for a spinning Kerr BH cannot be made conformally 
flat (at least for a broad class of foliations) lfl4l . and assuming CFA for spacetimes with 
rapidly-spinning BH has been found to introduce significant junk radiation lfT5l [Toll . The 
only way to accurately evolve general BH spacetimes is using full GR. Another limitation 
of Newtonian and CFA calculations is that gravitational waves must be estimated using 
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the quadrupole approximation, while in GR simulations the waves come directly from the 
spacetime evolution. 

The first challenge for a full GR evolution is to acquire valid initial data. This must 
both satisfy the ADM constraint equations and represent an astrophysically reasonable late- 
inspiral binary configuration. During the inspiral, the infall timescale is much longer than 
the orbital period. Also, gravitational radiation will damp away any orbital eccentricity 
before the separation becomes small. The exception might be mergers in dense stellar 
clusters, where 3-body interactions can be important, but these are expected to be a minority 
of NSNS and BHNS binaries ifTTl [T8l . (This expectation has recently been challenged, 
however lfT9ll .) Therefore, circular orbit is usually taken to be be a good first approximation 
for the late inspiral. Combining the circular orbit assumption with CFA, one can construct 
quasi-equilibrium profiles of the spacetime and the NS fluid. By computing these profiles 
for many binary separations, one can construct a sequence of "snapshots" that tell the history 
of the binary's inspiral. [To make each snapshot correspond to the same binary, the baryon 
mass of the NS(s) and irreducible mass of the BH are fixed.] During the late inspiral, the 
quasi-equilibruim approximation breaks down and full time evolutions are needed. A quasi- 
equlibrium snapshot of the late inspiral is then used as initial data for the evolution. Neglecting 
the initial radial (i.e. infall) velocity is known to produce spurious eccentricity in the orbit, but 
this can be removed by adding back an initial infall lt20ll2Tll22ll . CFA produces an initial burst 
of spurious "junk" gravitational radiation. Junk radiation can make it particularly difficult to 
create initial data with a rapidly spinning BH. This difficulty can be reduced by using a Kerr- 
Schild background metric, rather than CFA ll23ll2TT| . Some junk radiation will still be present, 
however, because the chosen conformal spatial metric still does not correspond to that of a 
true inspiraling binary. 

Evolving the metric forward in time requires choosing a formulation of GR and choosing 
a gauge. There are currently two formulations of GR used for NSNS/BHNS simulations: 
Baumgarte-Shapiro-Shibata-Nakamura (BSSN) (MUSI and generalized harmonic (GH) J26l 
l27l l28l |29l . In both formulations, the "trick" is to promote certain first spatial derivatives 
of the metric to the status of independent variables; this turns the evolution equations for 
metric components into wave-like equations. The gauge conditions fix the evolution of the 
coordinates. For BSSN, this is done by choosing the lapse a and shift vector (3 % evolution 
equations; the "1+log" lapse and "Gamma-driver" shift have proven immensely successful. 
With these gauges, BSSN codes can stably evolve spacetimes with BH singularities inside 
the computational domain ll30l l3"TTl . This treatment of BH spacetimes is called the "moving 
puncture" approach because BHs can freely move through the grid, a feature which has made 
it possible to do long-term, stable evolutions of binary black hole mergers with BSSN Il32ll33l . 
For GH, the gauge is specified through the gauge source functions H a = — y^TT ay g T . Suitable 
choices for H a are still under development (e.g. Il3"4ll35l0 . and the moving puncture approach 
has not yet been implemented in GH evolutions. Instead, BH singularities are excised in 
these simulations, i.e. a region inside the horizon is cut out of the grid and replaced with 
an inner boundary. Pretorius has successfully simulated binary black hole mergers using GH 
with excision |[36l[34]|. The Caltech-Cornell-CITA group has also successfully used GH to 
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simulate a variety of binary BH mergers using high-accuracy spectral techniques |[371[38l . 

The state of the neutron star fluid at each point is given by six numbers: the baryon 
density p, temperature T, 3-velocity v\ and electron fraction Y e . The fluid is evolved using 
six conservation equations: baryon number conservation (pw M ) ;M = 0, energy-momentum 
conservation T 7 ^^ = 0, and lepton number conservation (pY e u^)-^ = Sl, where the lepton 
source term Sl comes from neutrino emission or absorption. 

The stress tensor T^ v depends on the pressure P and specific enthalpy h; these are 
determined by the equation of state (EoS), which gives P and h as functions of p, T, and Y e . 
The actual EoS of NS matter is not known, so several guesses and approximations are used 
for NSNS/BHNS simulations. (For reviews of the NS EoS and its astrophysical implications, 
see ll39l |4Q| .) One is to treat the NS as a polytrope: P = np 1+1 / n , where the polytropic 
index n is a constant. This EoS is not temperature-dependent, so evolutions with this EoS 
will not allow the NS to heat. A more common EoS choice is to use the polytropic law to 
set the initial data for a cold NS and then evolve using the T-law P = (T — l)pe, where e is 
the specific internal energy and T = 1 + 1/n. The T-law is equivalent to the polytropic law 
as long as the fluid is continuous and isentropic, but it will allow heating when shocks are 
formed. Polytropic evolutions are still used to help gauge the importance of heating effects 
(e.g. ||4T|0 . A step towards more sophisticated EoS is to break up the pressure and internal 
energy into zero-temperature and thermal components, e.g. P = P co \d(p) + -Pthermai(p, T). 
Then the thermal part is set to a T-law [Pthermai = (r — l)pe t hermai], while the cold part 
can be a complicated function of p. Recently, piecewise-polytropic laws for the cold EoS 
have been proposed p2| and used ll43l . With only four parameters, this class of EoS can 
approximate any of the NS EoS proposed to date. It thus provides a good systematic way to 
cover the space of possible EoS. Other simulations use cold EoS based on nuclear physics 
calculations. Shibata and collaborators H4l I45ll22ll have evolved NSNS binaries in GR using 
the FPS EfO, SLy EZl, and APR flU cold EoS. The most realistic nuclear theory-based EoS 
have a general temperature and composition dependence. The tabulated Lattimer-Swesty H9l 
and Shen J50l E0 EoS have been used in Newtonian and CFA NSNS/BHNS simulations. 
Correct temperature dependence will probably be important for modeling mergers but not 
inspirals. Oechslin, Janka, and Marek ll52l have compared many of these EoS in the context 
of NSNS mergers. 

The transport of energy and lepton number by neutrinos can have important effects on the 
post-merger system. A few Newtonian NSNS/BHNS merger calculations have included these 
effects in approximate ways. The simplest way is through a leakage scheme ll53ll54l . These 
model neutrino cooling by removing energy from each gridpoint and adjusting Y e at a rate set 
to be some function of the fluid variables and the ^-optical depth r v . This function is chosen to 
reproduce rates given by local reaction timescales in //-transparent regions, while it reproduces 
rates given by the diffusion timescale in //-opaque regions. Note that the scheme is local: z/'s 
emitted in one part of the grid can't be absorbed in another part, a fundamental limitation 
of leakage schemes. Also, one must somehow estimate r v . A better but more complicated 
neutrino scheme, used in one study of the NSNS post-merger system ll55l is flux-limited 
diffusion ll56l I571l . In this case, the neutrino radiation intensity is evolved via a diffusion 
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equation, with fluxes limited so that free streaming is recovered in the //-transparent regime. 
Even these simulations do not solve the full multi-angle Boltzmann transport equation, a task 
beyond current numerical resources. Much less work has been done on //-transport in GR. De 
Villiers ll58l and the Illinois (UIUC) group ||59l have independently developed GR radiation 
transport schemes for the optically-thick limit, but these codes have not yet been applied to 
NSNS/BHNS binaries. 

Recently, some NSNS merger simulations have included magnetic fields |[60l[6Tl[62l[63ll . 
These simulations all work in the ideal MHD limit-a good approximation in most cases, given 
the high conductivities inside NSs. To investigate cases where resistive effects are important 
(in low-density, high-temperature regions), Palenzuela et al ||64| have developed, but not yet 
applied, a code to evolve the relativistic resistive MHD equations. 

Modeling the important physical effects is only one aspect of NSNS/BHNS simulations. 
There is also the computational challenge of adequately resolving the several relevant length 
scales of the problem: the gravitational wavelength, the stellar radius, and the length scales 
of various fluid and MHD instabilities. Some simulations (e.g. (651 [66l HH [68l H3 [52l [VOj) 
use Smoothed Particle Hydrodynamics (SPH). The representation of the fluid by particles can 
naturally adapt to changes in the fluid's size and shape. SPH has not yet been implemented in 
full GR. Most codes use grids to represent the fluid and metric. The problem of multiple scales 
is sometimes handled using Berger-Oliger moving-box-in-box adaptive mesh refinement 
(AMR). This technique is used by the Carpet module of the Cactus code 11711 . which is used 
in the UIUC EH and Whisky [|71 GRMHD codes. AMR is also used in the SACRA El 
code and by the BYU/LSU/LIU group's code [|75l . Testing is also being done on the use of 
multipatch ("cubed spheres") grids in place of Cartesian boxes 117611771178117911 . Multipatch 
grids allow one to extend the grid outer boundary R at a cost proportional to R rather than 
R 3 . Finally, the Caltech-Cornell-CITA (CCC) group uses a two-grid technique, inspired 
by the CoCoNuT stellar collapse code [80]. The CCC code SpEC evolves the metric 
pseudospectrally on one grid while evolving the fluid using standard finite-volume shock- 
capturing techniques on a second grid [18T1 . This code's main advantage is that the fluid grid 
needs only extend as far as the NS matter, rather than out to the gravitational wave zone. 

A final challenge in NSNS/BHNS modeling is that of covering the entire physically 
interesting range of binary parameters. NS masses may vary by tens of percents, while BH 
masses can vary widely. In addition, the BH can have a spin, with Kerr spin parameter o/Mbh 
anywhere from zero to one, oriented in any direction. The NS can also have a spin. However, 
the NS viscosity is thought to be too low for tidal effects to hold the NS in corotation If8~2ll83ll . 
so by the late inspiral, the orbital motion should dominate over the NS spin in most cases. 
Therefore, many simulations assume the NS to be initially irrotational. The effects of non- 
zero temperature (from tidal heating) lf8~2l [84l and magnetic fields |[62l [631 are expected to 
be negligible before the merger, and the pre-merger composition of the NS is fixed by (3- 
equilibrium. Therefore, the main physical parameters to be varied are the masses and (for 
BHNS binaries) the BH spin. Also, the EoS, since it is not known, must be treated as another 
parameter (or, rather, set of parameters) to be varied. 
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2. Neutron star-neutron star mergers 

2.1. History 

The first numerical simulations of NSNS mergers used Newtonian physics and polytropic 
EoS ll85l l86l l87l l65l l66l l88l . These studies found that the binary merges into a massive 
star rotating rapidly and differentially, but, being Newtonian, they could not check for 
the possibility that this object collapses to a BH. Simulations using CFA gravity and 
the zero-temperature Mayle-Wilson EoS ll89ll were carried out by Wilson, Mathews, and 
Marronetti ||90l[T2l . These pioneering computations surprisingly predicted that the NSs can 
collapse individually to BHs before merging. However, they used an EoS with fairly low NS 
maximum mass, and they were plagued by an error in one of their equations [|9T1 . so individual 
NS collapse is no longer considered likely ||92l . 

After this, the field bifurcated-some groups concentrated on improving microphysics 
while retaining Newtonian gravity, and some groups concentrated on improving the treatment 
of gravity while retaining simplified microphysics. In the first category would be the 
simulations of Ruffert, Janka, and collaborators ll53l l93l |54| and those of Rosswog and 
collaborators ll68l l54l l94l . These simulations used the finite-temperature EoS of Lattimer- 
Swesty or Shen, and they used leakage models of neutrino cooling. Meanwhile, binary 
polytropes were evolved in CFA gravity by Oechslin, Rosswog, and Thielemann |95| and 
by Faber, Grandclement, and Rasio ll9~6ll . For full GR models, an important step was the 
development of accurate quasi-equilibrium configurations to serve as initial data 11971 l98l 
l99l 1 1001 11011 . Finally, the first fully GR simulations of NSNS merger were performed by 
Shibata and collaborators H102[|103lll04L Since this time, efforts have been made to combine 
relativistic gravity with realistic microphysics. Oechslin, Janka, and Marek ll52l 1 10511 used 
realistic NS EoS in CFA simulations, while Shibata, Taniguchi, and Uryu |]44l l45l modeled 
mergers in full GR using realistic zero-temperature EoS (plus T-law thermal components). 
In the past few years, magnetic field evolution has been added to both Newtonian lf60l and 
GR ||6T1 l62l 1631 NSNS simulations. Also, the accuracy of these simulations has been greatly 
improved through the use of AMR H106ll4Tll74ll . Finally, efforts have been made to study the 
effects of MHD instabilities H1071 110811 and neutrino energy transport 111091 11101 11111 1551 on 
the post-merger system. The numerical studies to date can be combined to form a coherent 
picture, our current best guess, of what happens when two neutron stars collide. 

2.2. The emerging picture 

When two NS merge, they form a massive remnant. The first question to answer is "does 
the remnant collapse to a black hole?" It might seem like it should. There is a maximum 
mass M<rov,max for nonrotating NS, above which the star is dynamically unstable. Mrov.max 
depends on the unknown EoS; predictions of its value fall in the range 1.5-2.5M [1391 . Two 
1.4M stars will merge into a roughly 2.8M Q object. This is significantly above even high 
predictions of M T ov,max> so collapse seems likely. However, rotation and, to a lesser extent, 
shock heating increase the maximum mass. A star uniformly rotating at its mass-shedding 



Numerical relativity confronts compact neutron star binaries: a review and status report 7 

limit (above which matter would be ejected from the equator) can be stable with a mass 
up to M sup « 1.2M TO v,max> due to centrifugal support H112II . Rigidly rotating stars with 
Mrov,max < M < M sup are sometimes called supramassive. Greater centrifugal support, and 
hence greater maximum mass, can be achieved if the star rotates differentially, with the center 
having a significantly higher angular speed f2 than the equator H113[I114| . A NS supported 
above M sup by differential rotation is called a hypermassive neutron star (HMNS). Numerical 
experiments in GR have demonstrated that HMNS are stable on dynamical timescales H115I . 
NSNS remnants are found by simulations to have angular speeds several times higher at the 
center than near the equator H102H . so they are good HMNS candidates. 

GR simulations find that, after the merger, a NSNS remnant may either collapse promptly 
(i.e. on a dynamical timescale) to a BH, or it may survive as a HMNS. The determining factor 
is the mass of the NSNS system Mnsns- If Mnsns is above a threshold mass M t h, the remnant 
collapses promptly; if Mnsns < Mh, it forms a HMNS. All of the GR codes are in agreement 
on this basic point dSH E31 HI BD EH ■ For T = 2 polytropes, M th « 1.7M TOV ,max MM, 
but for more realistic, stiffer EoS, M th is found to be about 1.3 — 1.35M T ov,max B51 - For 
Mrov,max = 2.1M , the threshold is M th = 2.7 - 2.8Af Q . This is a typical NSNS binary 
mass, so it is quite possible that both prompt collapse and HMNS formation occur in nature. 

The prompt collapse case can be a viable GRB engine only if the post-collapse BH is 
surrounded by a massive accretion disk. For equal-mass binaries, the disk mass is usually 
very small 0.01M Q ) and the disk is very thin ll45l . making these systems poor candidates 
for producing GRBs. On the other hand, simulations find that the merger of unequal mass 
binaries proceeds much differently and can produce large disks ll65l 11161 11041 l45l |22| . In 
unequal mass binaries, the lower-mass NS can fill its Roche lobe before merger, so that it 
sheds matter both inward onto the more massive NS and outward into a tidal tail. The tail 
then falls back and contributes to the accretion disk around the BH. For a binary with mass 
ratio around 3:4, the disk can be 0.01M Q or larger, an excellent setup for a GRB ||4~5] . In 
addition to the mass ratio, the disk mass also depends on the mass of the system-binaries with 
M NSNS slightly above M th produce more massive disks than those produced by more massive 
binaries lfT04l l22ll. 

For M N sns close to M th , the remnant survives as a HMNS, but only for a few 
milliseconds. This case was mentioned by Shibata, Taniguchi, and Uryu [44] and studied 
in detail by Baiotti, Giacomazzo, and Rezzolla pT|. In the latter's simulation, the NSs merge, 
but then a bar-mode instability causes the two cores to split. The cores merge and split four 
times over a time of 8ms. During this time, the remnant is highly dynamical and asymmetric, 
and so it radiates gravitational waves. This radiation carries away angular momentum, so that 
the remnant loses centrifugal support and collapses, leaving a BH surrounded by a massive 
(0.07M Q ) disk. 

For lower-mass systems, the product of the merger is a HMNS. Such objects 
are dynamically stable, but they are vulnerable on longer timescales to processes (like 
gravitational radiation or magnetic fields) that sap the centrifugal support in the remnant's 
core. If the star is axisymmetric, its rotation will not cause it to emit gravitational waves. 
However, for rapidly rotating stars with stiff EoS, the minimum energy configuration will be 
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ellipsoidal rather than spheroidal H117H . (A famous example would be the Maclaurin spheroid 
and Jacobi ellipsoid configurations for rapidly rotating incompressible stars H118II .) Stiff stars 
tend to take ellipsoidal shape if (3 = T/\W\, the ratio of rotational kinetic to gravitational 
potential energy, exceeds a critical value. Realistic NS EoS tend to be stiff (T ps 2.75) at 
high densities, and numerical simulations with Newtonian ll65l I119H . post-Newtonian HI 201 
and full GR flM BD 321 physics have indeed found that NSNS remnants with these EoS 
usually form ellipsoids. These bar-like deviations from axisymmetry are stronger for stiff er 
assumed EoS ll52~l . Such stars emit strong gravitational waves and lose angular momentum 
J. This might have the effect of causing the star to contract or of lowering the ellipticity, but 
simulations indicate that the former effect dominates, and the remnant remains ellipsoidal as 
it contracts. Eventually, the remnant will reach a critical state and collapse to a BH. From 
the J-loss rate, the time for this to happen is inferred to be 30-100ms. Such long timescales 
are difficult to simulate in 3D with current resources, so the first simulations stopped well 
before the collapse. Oechslin and Janka H105H and Shibata and Taniguchi ll45l estimated the 
mass of the disk around the post-collapse black hole based on the J distribution in the HMNS. 
Matter with specific angular momentum high enough for orbit outside the ISCO of the (not 
yet formed) BH was assumed to form the accretion disk. These calculations give disk masses 
of ~ O.IMq, high enough for GRBs. However, these estimates necessarily could not account 
for changes in the HMNS's J distribution prior to collapse. Finally, Baiotti, Giacomazzo, and 
Rezzolla ([4l| carried out long-time (~ 30ms) simulations using a cold polytrope EoS that, for 
the first time, evolved the NSNS remnant to the point of delayed collapse (in their case, 16ms 
after merger). They found a final disk mass of ~ O.O8M , quite close to the earlier estimates. 

Prompt collapse or ellipsoidal HMNS formation occur in most simulations, but other 
evolutionary paths are possible. If (3 of the remnant is below the critical value, the HMNS 
may form a spheroid which will eventually be driven to collapse by magnetic-driven, rather 
than gravitational wave-driven, processes (see below). This might happen for M NSNS near but 
below M t h for some EoS ll45l . Also, if Mnsns < M sup , the remnant might not collapse at all. 

Shocks during the merger heat the remnant to temperatures of order 10 MeV. At these 
densities and temperatures, the remnant will radiate primarily in neutrinos. The HMNS has 
//-optical depth of 10 2 — 10 4 , depending on the neutrino energy, and it will cool from radiation 
on a timescale of order 1 sec [|55l . The neutrino emission from NSNS remnants has been 
studied using leakage schemes by Ruffert, Janka, and collaborators ll53l l93l I121R and by 
Rosswog and Liebendorfer ll54ll . More recently, Dessart et al ll55l have modeled neutrino 
transport in the HMNS under the assumption of axisymmetry using the multi-group flux- 
limited-diffusion code VULCAN/2D. These three codes are able to extract the neutrino signal 
from NSNS coalescence. They all find that the v e emission dominates over the v e emission, 
an unsurprising result given the low Y e of the remnant. Dessart et al also find high v^lv T 
luminosities, annihilation releases some 10 49 — 10 50 erg s" 1 outside the remnant. However, 
neutrino heating also drives a wind from the HMNS surface. This wind dumps ~ 10~ 3 M Q s _1 
of matter into the polar regions, enough to baryon-load the poles and prevent a GRB from 
being produced while the HMNS remains ll55l . Dessart et al speculate that a similar wind 
may be produced by the disk after the HMNS collapses, but in this case a centrifugal barrier 



Numerical relativity confronts compact neutron star binaries: a review and status report 9 

might keep the polar regions baryon-clean, allowing a GRB. The //-driven wind also carries 
away some 10~ 3 M Q of neutron-rich material out of the HMNS system. GR may modify the 
above results, given that GR simulations tend to produce less massive tori around the HMNS 
than their Newtonian equivalents, and these tori contribute disproportionately to the neutrino 
emission ||54~1 . 

Newtonian lfT22l[T211 and CFA J52l simulations predict that 10~ 3 - 1O~ 2 M , depending 
on the EoS, of nuclear matter is ejected during NSNS mergers. The CFA simulations found 
two components to the ejecta: cold matter coming off the tidal tail and hot matter ejected from 
the surface of contact between the stars. As the expelled matter decompresses, neutrons are 
captured onto the heavy nuclei. The resulting nucleosynthesis has been modeled, with various 
approximations concerning the expansion rate and the initial composition and temperature, 
using r-process network calculations H1231 1124| . For suitable parameters, the calculated 
abundances agree quite well with solar abundance pattern for mass numbers A > 140, 
especially regarding the A = 195 peak. However, reliable estimates of the ejecta mass and 
composition require GR simulations with realistic EoS. 

The effect of the magnetic field in a NSNS binary is probably negligible prior to the 
merger for realistic field strengths ll62l l63l . For very large fields (> 10 16 G), magnetic 
tension can suppress tidal deformation of the NSs and delay the merger lloTl l63l . When 
the stars touch, a shear layer is formed at their interface. The shear layer is Kelvin-Helmholtz 
unstable, causing the flow to curl up into vortex rolls. Inside these vortex rolls, the magnetic 
field can be quickly amplified by orders of magnitude. This process was first modeled using 
Newtonian MHD by Price and Rosswog ll60l . From a 10 12 G pre-merger field in each NS, 
they found that the field in vortex rolls reached 10 15 G in 1ms. They speculate that the field 
will reach equipartition in small vortices. These high-field pockets of matter will become 
buoyant, float up, and produce relativistic blasts when they hit the NS surface. If so, these 
blasts might help provide the energy for a GRB, and magnetic pressure might help deflect the 
neutrino-driven baryon wind. Kelvin-Helmholtz-driven field amplification has also been seen 
in GR simulations lloTl l63l . In the simulations of Giacomazzo et al ll63l . the magnetic field 
saturates at more modest values (e.g. a 10 12 G field grows by less than an order of magnitude). 
The discrepancy with the result of Price and Rosswog reflects the difficulty in numerically 
capturing the smaller- scale effects introduced by the Kelvin-Helmholtz instability. 

Because the angular speed of the HMNS decreases with radius, it is unstable to the 
magnetorotational instability (MRI) H125L The MRI produces small-scale turbulence that 
redistributes angular momentum outward, causing the core to slow, contract, and eventually 
collapse. This process has been modeled in GR MHD under the assumption of axisymmetry 
by a collaboration of the UIUC group and Shibata [fT07l[T08l[T26ll . They found that the HMNS 
collapse leaves a BH surrounded by a massive accretion disk and baryon-clear poles. MRI- 
induced delayed collapse is likely the fate of spheroidal HMNS and ellipsoidal HMNS with 
very strong magnetic field. 

Much of the interest in NSNS binaries is motivated by the prospect of extracting 
information about NS structure from inspiral and merger gravitational waveforms. The 
inspiral waveform may be particularly important for these purposes, because it includes 
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the frequencies at which Advanced LIGO will be most sensitive. The early, low-frequency 
inspiral waveform can be used to constrain the NS tidal Love number H127H . Faber et al H128II 
point out that the energy spectrum of the late inspiral waveform can constrain the NS radius. 
Read et al II431 have looked for EoS signatures in the late inspiral waveform by performing 
GR NSNS simulations and systematically varying the EoS. They conclude that the NS radius 
can be determined from the waveform to an accuracy of ~1 km for an event at 100 Mpc. 

The merger, HMNS, and BH ringdown waveforms all have frequencies above 1kHz and 
so are more difficult to detect, but the mass and EoS-related differences become quite dramatic 
in them. If Af NSNS > M th , the inspiral signal is followed by a merger waveform at frequecies 
1kHz < / < 3kHz, which in turn is followed by a (probably undetectable) BH ringdown 
signal at 6.5-7kHz ||22l . If M NSNS > M th , the spectrum has another peak at ~3kHz due to 
the radiation from the ellipsoidal remnant. A HMNS lasting ~ 50ms located ~ 50Mpc away 
could be detected by Advanced LIGO with a signal-to-noise ratio of around 3 (1 12911 . On its 
own, this would be a weak signal, but it will be preceeded by a more easily-detectable inspiral 
signal. The EoS strongly affects the compaction of the HMNS, and therefore its moment of 
inertia and its rotation frequency. Therefore, the frequency of the post-merger signal peak 
contains information on the EoS, particularly if Mnsns can be extracted from the inspiral 
signal lfT29irT30H . 

3. Black hole-neutron star mergers 

3.1. History 

The first simulations of BHNS mergers were carried out a decade ago by Lee and 
Kluzniak ll67l 1 1 3 1 L 11321 1133H . They used Newtonian physics, with the NS modeled as a 
polytrope and the BH modeled as a point mass. Next, simulations in Newtonian physics with 
nuclear- theory EoS were performed IU34U1351 . and the Newtonian potential was replaced by 
a Paczyhsky-Wiita potential H136U137l l. Faber et al ll69ll took a step toward GR by simulating 
extreme mass-ratio binaries using CFA physics. Loffler, Rezzolla, and Ansorg H138H modeled 
a head-on collision between a black hole and a neutron star in full GR. GR simulations 
of orbiting binaries became possible after quasi-equilibrium initial data was successfully 
generated ffHH EES [HTj \TM EES [2TJ- The first full-GR simulations of BHNS mergers 
were carried out by Shibata and Uryu |)144l I145H . followed quickly by the UIUC H146II 
and CCC ifSTI groups. These simulations all used polytropic EoS and only considered the 
case of zero initial BH spin. Meanwhile, Rantsiou et al [|70l modeled BHNS mergers with 
spinning BH using a Kerr spacetime and an approximate treatment for the NS self-gravity. 
Their findings suggesting the importance of BH spin were confirmed when the UIUC group 
simulated the merger of several BHNS systems with nonzero initial BH spin in full GR ll72l . 

3.2. The emerging picture 

The mass of the black hole, M B h, can vary widely from one BHNS system to another. For 
very high M B h, the NS behaves like a test particle: it inspirals slowly until it reaches the ISCO 
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separation g? isco , at which point it plunges into the BH and is swallowed whole. For M BH of 
only a few times the neutron star mass M NS , the NS will overflow its Roche lobe at a radius 
efaisr before a plunge takes place. In this case, which resembles the unequal-mass NSNS case, 
the NS is tidally disrupted while outside the BH. Matter flows into the hole and outward into 
an extended tidal tail. Thus, we expect tidal disruption if cfaisr > g?isco> an d a plunge into 
the BH otherwise. In addition to the mass ratio, ofaisr depends on the NS radius -Rns- A more 
compact star can survive stronger tidal forces and will disrupt closer to the BH. For expected 
i?NS> the critical mass ratio is about 4:1 [|142] . For higher M BH or smaller i? NS , one would 
expect the NS to be swallowed whole; for lower M BH or larger i? NS , tidal disruption and mass 
transfer are expected. 

The latest GR simulations by Shibata et al H 14711 confirm these expectations. Using a 
T = 2 polytropic EoS and fixing the initial BH spin to zero, they vary the mass ratio over the 
range 1.5:1 to 5:1, and they vary i?NS over the range ll-13km. For mass ratio 5:1, the NS is 
swallowed whole, leaving no disk. The gravitational waveform contains inspiral, merger, and 
BH ringdown signals which are similar to a binary BH merger of the same masses. The final 
kick velocity is also similar to the binary BH case. For mass ratios below 3:1, the NS disrupts 
before falling into the hole. The disruption spreads out the NS matter, so that the matter 
is more symmetrically distributed around the BH as it accretes, and the gravitational wave 
emission of the system greatly diminishes. This manifests itself in the Fourier spectrum of 
the waveform as a steep decline in amplitude above a cutoff frequency / cut ~ 1.3 /tidal, where 
/tidal is the wave frequency at d disi . Since d disr depends on -Rns, /cut contains EoS information. 
If the star disrupts during its plunge phase, the waveform will have inspiral and merger signals, 
but the ringdown will be significantly weakened. If the the NS disrupts during the inspiral, 
then the merger part of the waveform is also suppressed. Disruption of the NS outside the BH 
also results in much lower kick velocities for the final hole than occur in the equivalent binary 
BH cases. Shibata et al H147L Etienne et al (114611 . and Duez et al ||8T1 all find that, under 
the assumption of T = 2 EoS and nonspinning BH, BHNS mergers usually produce small 
accretion disks. Disk masses of ~ 1O~ 2 M were obtained by Shibata et al 1114711 for very low 
mass ratios (e.g. 2:1) and large NS radii 14km) 1114711 : in all other cases, the disk masses 
were much lower. In their most recent paper 11721 . the UIUC group find somewhat higher 
disk masses: they find a 0.06M o disk for a case with a 3:1 mass ratio, -Rns ~ 14km, and no 
initial BH spin, a case that had not produced a significant disk in previous studies II 14611 147ll . 
Sources of inaccuracy in these studies include limited resolution and inadequate treatment of 
low-density regions. (See [1721 for a discussion of the latter issue.) 

Things are different when the BH is spinning, because <iisco f° r prograde orbits around 
a spinning Kerr BH is smaller than c^sco f° r a nonspinning hole of the same mass. Therefore, 
a NS would be expected to get closer to the BH and experience greater tidal forces before 
plunging. Etienne et al 11721 have recently simulated several BHNS systems with nonzero 
BH spin aligned with the orbital rotation axis. For a mass ratio of 3:1 and a T = 2 EoS, 
they consider an aligned spin case with a/M B H = 0.75 and an anti-aligned spin case with 
a/M B H = 0.5. They find that the binary with aligned BH spin merges to produce a disk with 
mass O.2M and temperature ~ 4MeV, a promising setup for a GRB. Aligned BH spin may 
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also make tidal disruption possible in cases where the BH is much more massive than the NS. 
An affine quasi-equilibrium model by Ferrari, Gualtieri, and Pannarale H148I predicts that, for 
a BH spin of o/Mbh = 0.9, NS will be disrupted for mass ratios up to 10: 1 to 20: 1, depending 
on the NS mass and EoS. 

The merger can be strongly affected by the NS EoS. The EoS determines -Rns> the effects 
of which have already been discussed. The EoS also determines how i? NS changes with the 
NS mass, which can affect the character of the mass transfer. If a small transfer of mass 
leads the star to expand relative to its Roche lobe, mass transfer will be unstable. If the star 
contracts into its Roche lobe, mass transfer might be stable or turn on and off. In general, 
mass transfer tends to be more stable for stiffer EoS (higher Y) because this corresponds to a 
greater tendency for a star to shrink when it loses mass [1 14911 . 

Attempts have been made to estimate the effects of NS EoS in the context of BHNS 
simulations using Newtonian gravity. Lee and Kluzniak ll67l 1 1 3 1 H and Lee H1321 I133B 
considered poly tropes with Y between 5/3 and 3. Janka, et al HI 3411 used the Lattimer- 
Swesty EoS BUI, while Rosswog, Speith, and Wynn lfT35ll used the Shen EoS ||50l |5H. 
These simulations showed large qualitative differences for different EoS assumptions. For 
Lattimer-Swesty nuclear matter, the NS disrupts in one mass transfer event, and a large post- 
merger disk is created. For Shen nuclear matter, a NS core can survive multiple mass transfer 
events, and the postmerger disk is much smaller. There are indications that the mass transfer 
is much less stable in GR, and the differences between the cases of stiff and soft EoS is not 
nearly so dramatic. The use of GR-mimicking potentials [TTQl ITTIl tends to eliminate episodic 
mass transfer HI 361 1 1371 . Simulations of large mass-ratio cases using CFA also find generally 
unstable mass transfer ||69l . Most recently, the CCC group has evolved BHNS binaries with 
3:1 mass ratio, -Rns = 14km, and a/M B H = 0.5 for a variety of EoS, including polytropic 
(T = 2 and Y = 2.75) and Shen H150II . They find that the NS is disrupted in one mass transfer 
event in every case, and the disk masses are all comparable (0.1-0.2M o ), although stiffer EoS 
do produce larger, longer- lived tidal tails. 

When tidal disruption decompresses the NS matter, neutrons and protons combine to 
form heavy nuclei, which heats and thickens the tidal tail H1351I1501 . Metzger et al H151II point 
out that, if matter from the tail remains in marginally bound orbit for around one second, even 
more energy may be released by r-process nucleosynthesis, perhaps enough to unbind this 
matter. Assuming a hot disk is sometimes formed, as in the simulations of Janka et al H134I . 
the nucleosynthesis in an outflow from such a disk has been calculated by Surman et al H152L 
They find that r-process nucleosynthesis does occur in such outflows. The likely BHNS 
merger rate and ejecta mass make it unlikely that this is the main source of these elements, 
however. 

4. The post-merger black hole plus disk system 

It appears that a massive accretion disk (M disk ~ 10~ 2 — 1O~ 1 M ) can be produced by 
either a NSNS or a BHNS merger. Such disks are initially dense (p ~ 10 10 — 10 12 g cm -3 ), 
hot (T ~ 1 — lOMeV), and neutron-rich (Y e ~ 0.1). Angular momentum transport, most 
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likely driven by MRI-induced turbulence II 12511 , causes matter to accrete into the black hole at 
enormous "hyperaccretion" rates of 0.1 M s -1 to lOM s -1 . 

Steady- state models of hyperaccretion onto black holes were constructed by Popham, 
Woosley, and Fryer II 15311 . Their models were one-dimensional (axisymmetric and vertically 
integrated) but incorporated GR using a Kerr background spacetime. Following the Shakura- 
Sunyaev prescription 1115411 . the angular momentum transport was modeled by adding a 
shear viscosity of strength v = olc^/VLk, where c s is the sound speed, VL K is the Keplerian 
angular velocity, and a is a (unknown) dimensionless parameter. Because neutrino emission 
was found to be the dominant cooling process, such disks are sometimes called "neutrino 
dominated advection flows" (NDAFs). Other steady-state, ct-viscosity NDAF models soon 
followed H1551 11561 11571 . with the most important improvement being the inclusion of 
neutrino opacity effects H158II . 

Time evolutions of NDAFs under the influence of a-viscosity have been carried in 
ID lfT59ll . 2D (axisymmetry) [fT60l [T6T1 fTTQll and 3D ifTOfTTTTl . The 2D and 3D simulations 
used as their initial conditions the final state of numerical BHNS merger simulations, and they 
treated neutrino processes via leakage schemes. Shibata, Sekiguchi, and Takahashi HI 621 have 
taken an important step by including magnetic fields in their NDAF evolutions. They evolve 
the MHD equations in GR but assume axisymmetry. This allows their code to capture effects 
related to the MRI, eliminating the need for an a viscosity. 

The stresses that drive angular momentum transport also cause energy to be dissipated 
as heat, raising the temperature of the disk to T ~ lOMeV. Photons in the NDAF are trapped, 
so the disk cannot cool by photon radiation. In the inner ~ 10 3 km, the disk is dense and 
hot enough to cool by neutrino emission. (According to most merger simulations, this would 
include the entire disk.) In sufficiently massive disks, for which the density reaches p > 10 11 g 
cm -3 , neutrinos can become trapped in the inner region, and thermal energy is advected into 
the hole rather than radiated. For M > lM Q s _1 , neutrino trapping significantly limits the 
radiated energy that might power a GRB H158II . 

To power a GRB, some 10 50 erg s _1 must be transferred from the BH-disk system to 
a low-mass, ultrarelativistic outflow. NDAF simulations have investigated the possibility 
that this energy is provided by neutrino emission. They find that the neutrino luminosity 
L v increases steeply with disk mass M d i S k, viscous strength a, and corotating BH spin 
a/M BH Ifl53l mTl [1571 fl62ll . In favorable circumstances (M disk ~ 0.1 M & , a ~ 0.1, 
a/M BH > 0.5), L u ~ 10 53 erg s _1 can be released, of which a few percent could be converted 
to a pair-photon fireball by vv annihilation ||111| . This is sufficient to explain observed short 
GRB energies if the outflow is collimated into a ~1% angle. GRMHD simulations of NDAFs 
also find L v ~ 10 53 erg s^ 1 , indicating that the "true" effective a is around 0.01-0.1 H162H . 
These simulations also found strong variations in L v on millisecond timescales, in agreement 
with the observed variability of GRBs. 
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5. Future work 

The needed future improvements in NSNS/BHNS merger modeling fall into three categories: 
improved microphysics, improved coverage of parameter space, and improved numerical 
accuracy. For GRB modeling, accurate microphysics is the most pressing need. Numerical 
simulations aiming to determine whether or not a sufficiently energetic, sufficiently 
ultrarelativistic outflow is generated will probably need GR, MHD, and neutrino cooling, 
heating, and scattering. A full solution to the neutrino transport problem is currently not 
feasible, but we may hope that approximate treatments will capture the essential physics, such 
as vv energy deposition and //-driven winds. For the other main goal of providing gravitational 
waveform catalogs for LIGO data analysis, covering the NSNS and BHNS parameter space 
is the most important issue. The ongoing studies of the effect of the equation of state on the 
waveforms using parameterized EoS are obviously vital. Another poorly studied variable is 
the BH spin in BHNS binaries. So far, only a few cases have been studied in GR, and all 
of them belonged to the special case of aligned spin. Numerically, one major challenge is to 
improve resolution to so that turbulent flows are adequately treated. Another is to extend the 
evolution time so that more NSNS mergers can be evolved to the point of delayed collapse. 
Fortunately, each of these issues is being pursued by at least one research group. 
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